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GENERAL  DISCUSSION 


The  FAST2D  code,  a  two-dimensional  hydrocode  developed 
by  the  Laboratory  for  Computational  Physics  at  NRL,  has  been 
modified  and  used  to  model  homogeneous  flows.  The  speed  of 
the  transport  algorithm  which  advances  in  time  the  density, 
momentum,  and  energy  has  been  improved  by  a  factor  of  two. 

This  improved  code  has  been  used  to  model  the  transition  to 
turbulence  for  shear  flow. 

When  fluids  of  different  densities  flow  at  different 
relative  velocities,  one  stream  above  another,  an  instability 
develops.  The  vorticity  of  a  slab  with  a  slight  initial  dis¬ 
turbance  is  positive  for  a  positive  velocity  leading  to  a 
migration  of  positive  vorticity.  Where  vorticity  is  displaced 
downwards  a  positive  x  component  is  induced.  This  accumulation 
will  induce  clockwise  velocities  which  amplify  the  initial  dis¬ 
turbance  leading  to  the  instability. 

To  model  such  a  flow  system  a  computational  grid  was 
chosen  200x64  in  x  and  y,  respectively.  Zone  sizes  of 
2.8cm  in  x  and  1.0cm  in  y  were  used.  Downstream  the  grid 
was  stretched  to  obtain  maximum  length.  To  insure  against 
spurious  signals  propagating  upstream  a  special  boundary 
condition  was  implemented.  The  momentum  outflow  was  compared 
to  the  momentum  inflow  and  adjusted  to  obtain  a  balanced  flow. 
This  boundary  condition  was  checked  for  several  flow  velocities 
to  insure  a  consistent  answer. 


The  outflow  boundary  value  for  the  pressure  was  chosen 
in  the  following  manner.  To  insure  that  the  average  pressure 
does  not  greatly  deviate  from  the  ambient  pressure,  the  out¬ 
flow  pressure  must  act  so  as  to  restore  the  computed  pressure 
to  ambient.  Therefore,  we  assert  that  the  condition 

P  =  P 

outflow  ambient 


would  be  an  overly  harsh  restoring  condition,  if,  for  exam- 
pie  P  »  1.05  Pamblent  or  P  <  0.95  Pambient.  The  gradient 


P-P 


ambient 


/AY  has  associated  with  it  some  characteristic 


scale  length.  If  l  is  the  width  of  our  system  such  a  length 
might  be  1/2,  then 


outflow 


¥  (pambient-p1  *  P 


Take  the  case  where  P  >  Pam5ient,  then  Poutflow  <  P  which 
creates  a  positive  gradient  allowing  flow  out.  The  case 
where  P  <  Pambient  allows  PQutflow  >  P  causing  a  negative 


An  important  part  of  studying  the  development  of  the 
instability  is  analyzing  the  data.  The  graphics  routines 


were  modified  and  enhanced  to  establish  a  color  presentation 
capability.  A  shadow  plotter  was  used  as  the  basis  for  the 
density  plotter.  A  vector  velocity  plotter  was  modified  to 
present  arrows  in  color.  Figure  1  shows  the  resulting 
display  of  density  and  velocity  for  the  problem  at  1.72 
milliseconds . 

Several  preliminary  calculations  were  completed  to  test 
the  overall  production  ability  of  the  code  with  the  new 
boundary  conditions  and  the  new  data  presentation.  Study 
of  these  co-flowing  problems  requires  many  sets  of  con¬ 
figurations,  and  thus  the  ability  to  produce  a  complete 
simulation  every  five  to  six  days  including  color  processing 
time.  This  requirement  has  been  met  and  the  code  is  now- 
ready  for  that  use. 

For  the  subsonic  reactive  flow  part  of  the  program, 
the  following  three  tasks  were  performed: 

(1)  the  one-dimensional  time  -  dependent  flame  model 
was  adapted  to  include  a  capability  for  modeling  the 
thermal  diffusion  phenomena, 

(2)  a  study  of  the  effects  of  thermal  diffusion  on 
the  burning  velocity  of  hydrogen-air  mixtures  was 
carried  out ,  and 

(3)  the  flame  code  has  been  documented  as  NRL 
Memorandum  Report  4910. 


Subsonic  Reactive  Flows 


In  the  past'*',  SAI  in  collaboration  with  NRL  had 
modified  and  used  the  NRL  Flame  code  to  study  the  ignition 
and  quenching  of  pre-mixed  gases0  and  for  the  determination 

4 

of  the  burning  velocity  of  fuel-air  mixtures  .  It  was  also 
4 

suggested  earlier  that  the  thermal  diffusion  phenomena 
might  play  a  significant  role  in  the  burning  of  hydrogen- 
air  mixtures.  During  this  contract  year,  the  flame  code  has 
been  extended  to  include  the  thermal  diffusion  phenomena  and 
the  effect  of  thermal  diffusion  on  the  burning  velocity  of 
a  stoichiometric  hydrogen-air  mixture  has  also  been  studied. 
A  detailed  account  of  this  work  is  included  in  this  report 
in  Appendix  A  and  a  brief  description  is  given  below. 

Thermal  diffusion  contributes  an  additional  source 
term  [see  Eq.  (8)  in  App .  A]  to  the  diffusion  equations 
[Eq.  (7)  in  App.  A]  and  in  this  way  changes  the  diffusion 
velocity  of  the  various  species  in  the  system.  It  also 
modifies  the  energy  equation  [Eq.  (4)  in  App.  A]  by  con¬ 
tributing  a  term  to  the  heat  flux  vector  [Eq.  (5)  in  App. A]. 
The  evaluation  of  these  terms  and  the  solution  of  the  dif¬ 
fusion  and  energy  equations  is  discussed  in  detail  in 
Appendix  A.  In  the  discussion  below,  the  source  terms  in 
the  diffusion  equations  due  to  thermal  diffusion  will  be 
represented  as  and  those  due  to  ordinary  diffusion 

(due  to  concentration  gradients)  as  Sgp. 
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In  order  to  evaluate  the  importance  of  thermal  diffusion, 
the  burning  velocity  of  a  hydrogen-air  mixture  was  determined 
with  and  without  considering  the  source  terms  due  to  thermal 
diffusion.  The  net  effect  of  thermal  diffusion  is  to  lower 
the  burning  velocity  of  a  stoichiometric  hydrogen-air  mixture 
by  about  1%.  In  order  to  understand  the  reasons  for  this 
effect,  the  spatial  variation  of  the  terms  SqD  and  for 

the  various  species  have  been  calculated  for  a  steady  propa¬ 
gating  flame.  For  the  species  09  ,  the  term  is  opposite 

in  sign  to  Sq^.  A  negative  value  for  implies  that  the 

oxygen  molecules  are  diffusing  from  the  hotter  region  to 
the  cooler  region.  Since  the  sum  of  the  two  source  terms 
determines  the  diffusion  velocity,  the  effect  of  thermal 
diffusion  is  to  slow  the  diffusion  of  0-,  molecules  into 
the  higher  temperature  region  of  the  flame.  The  effect  of 
this  on  the  burning  velocity  is  not  obvious  since  it  depends 
on  the  details  of  the  reaction  mechanism  and  the  distri¬ 
bution  of  other  species.  For  the  species  N9 ,  the  term  S^j-, 
is  comparable  to  the  term  Sq^  and  the  effect  of  thermal 
diffusion  is  to  enhance  the  movement  of  N?  from  the  hotter 
region  of  the  flame  to  the  cooler  region.  For  the  species 
H?,  the  effect  of  thermal  diffusion  is  similar  to  that  of 
ordinary  diffusion  and  it  enhances  the  movement  of  H?  from 
the  cooler  region  to  the  hotter  region. 
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The  term  Sq^  for  the  species  H  is  particularly 
interesting.  It  is  observed  that  the  H  atoms  in  a  certain 
region  of  the  flame  (the  region  with  T  >  1400  K  for  the  flame 
studied  here)  have  a  tendency  to  diffuse  into  the  hotter 
region  while  the  rest  (in  the  region  with  T  <  1400  K)  diffuse 
into  the  cooler  region.  The  diffusion  of  H  atoms  into  the 
cooler  region  is  very  important  in  the  propagation  of  the 
flame.  The  effect  of  thermal  diffusion  is  to  delay  and 
decrease  the  diffusion  of  H  atoms  into  the  cooler  region. 

Since  the  H  atoms  are  produced  primarily  in  the  high  temper¬ 
ature  region,  the  thermal  diffusion  of  H  atoms  can  signifi¬ 
cantly  retard  the  rate  of  propagation  of  the  flume.  This 
may  be  the  main  reason  for  the  observed  reduction  in  the 
burning  velocity  when  the  effects  of  thermal  diffusion  were 
included.  This  has  been  discussed  in  more  detail  in 
Section  5.2  of  Appendix  A. 

The  model  used  for  the  flame  studies  discussed  above 
and  elsewhere^  ^  is  a  one-dimensional,  time-dependent, 
Lagrangian  numerical  model.  The  model  incorporates  a 
number  of  new  approaches  and  algorithms  which  have  been 
tested  by  comparisons  to  less  complex  or  analytic  solutions 
and  by  comparisons  to  experimental  data.  These  new'  approaches 
and  algorithms  along  with  the  input  parameters  required  by  the 
model  have  been  discussed  in  the  NRL  Memorandum  Report  4910 


b 


entitled  '  A  One -Dimensional  Time -Dependent  >*ode]  for  Flame 
Initiation,  Propagation  and  Quenching."  This  memorandum 
report  has  been  included  in  this  report  as  Appendix  A. 
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A  ONE*DIMENSIONAL  TIME-DEPENDENT  MODEL 
FOR  FLAME  INITIATION,  PROPAGATION  AND  QUENCHING 

1.  INTRODUCTION 

This  report,  describes  a  one-dimensional,  time-dependent ,  lagrangian 
numerical  model  developed  to  study  the  initiation,  propagation  and  quenching 
of  laminar  flames.  The  model  incorporates  a  number  of  nev  approaches  and 
algorithms  which  have  now  been  tested  by  comparisons  to  less  complex  or 
analytic  solutions  and  by  comparisons  to  experimental  data.  These  nev  ele¬ 
ments  include:  ATTNC  [  1 ]  an  implicit,  Lagrangian  method  for  solving  the 
convective  parts  of  the  conservation  equations;  DFL’JX  [2,3]  a  variable 
accuracy  algorithm  for  determining  diffusion  fluxes  without  having  to  invert 
matrices;  VS AIM ,  a  vectorized  version  of  the  ordinary  differential  equation 
solver,  CHEMSQ  ] 1 , 5 ] ;  and  a  new  method  for  treating  an  open  boundary  in  an 
implicit,  Lagrangian  calculation.  An  asymptotic  coupling  method  is  used  in 
conjunction  with  timestep  splitting  to  couple  the  various  physical  and  chemi¬ 
cal  processes.  This  approach  allows  the  use  of  entirely  different  algorithms 
for  the  physical  processes  represented  by  the  different  terms  in  the  conser¬ 
vation  ecuatior.s. 


Manuscript  submitted  July  26,  1982. 


*•  •-  •*  ' 


•vr  - 


P 


om  3AC\G~.C  mCS  ~ 


lasinar  52.2.2s  is  ihs  sarliss*  combustion  rroclsz  -0  ot 

studied  theoretically  vhich  required  the  simultaneous  consideration  of  bcth 
fluid  dynamics  a^d  chemical  kinetics  for  its  solution  [l“[.  The  problem  cf 
determining  the  propagation  velocity  cf  a  deflagatior.  wave  was  first  studied 
'ey  Mallard  and  leChatelier  f  11] .  Since  then  there  have  been  many  analytical 
attempts,  fe.g.  12,  13]  vitfc  varying  degrees  of  complexity  and  approximation , 
to  study  s imnle  steady-  flames.  Some  of  these  have  beer,  summarised  by 
Williams  !  10 ;  .  'whereas  these  analytical  techniques  are  extremely-  informa¬ 
tive,  they  cannot  satisfactorily  describe  the  detailed  structure  cf  flames 
since  they  all  include  very  approximate  chemical  kinetic  schemes  and  trans¬ 
port  properties.  For  this  one  has  to  resort  to  numerical  methods. 

One  of  the  first  attempts  to  solve  the  laminar  flame  problem  numerically 
vas  that  of  Mrschf elder  and  co-workers  [ihj.  They  formulated  the  unsteady 
flame  problem  as  a  system  of  three-dimensional  nonlinear  partial  differential 
equations  and  they  solved  the  one-dimensional  steady  flame  as  a  two  point 
boundary  value  problem.  They  used  approximate  methods  to  estimate  the  mass 
flow  rate  which  is  the  eigenvalue  of  the  two-point  boundary  value  problem  and 
then  used  a  numerical  shooting  method  to  obtain  the  temperature  and  species 
profiles.  Although  the  first  problem  they  studied  involved  only  single  step 
kinetics,  they  later  applied  the  same  solution  procedure  to  study  flames 
for  which  the  kinetics  involved  chain  reactions  !l5»l6]. 

In  l?5c,  Spalding  [IT]  applied  a  time -dependent  system  cf  nonlinear 


partial  differential  equations  using  explicit  finite-difference  techniques  no 
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methane,  oxygen,  nitrogen  fixtures  using  very  detailed  chemical  kinetic 
mechanisms  and  transport  properties. 

7r.e  kinetics  and  propagation  of  methane-air  flames  vas  studied  'ey  Erect 
et  al.  1 2 6 !  using  a  mixed  explicit-implicit  finite  difference  technique.  7r. 


diffusive  transport  t 


were  solves  explicit_y  ano  tne  Kinetic  terms 


solved  using  linearised  implicit  techniques.  -ledgian  ;2” ’  employed  a  meth: 
of  lines  technique  vhere  the  nonlinear  parabolic  initial-boundary  value 
problem  is  reduced  to  a  set  of  nonlinear  first  order  initial  value  problems. 
Sophisticated  initial  value  problem  integrators  were  used  to  solve  the 
resulting  ordinary  differential  equations  system.  Instead  of  using  finite 


approximations  to  discretize  the  spatial  derivative,  Margclis  .2 


used  a  spline  collocation  procedure  in  his  method  of  lines  approach, 
■'•estbrook  and  Dryer  [2?'  have  deduced  a  comprehensive  reaction  mechanism  fer 
methanol  oxidation  involving  eighty  four  reactions  using  an  implicit  finite- 
difference  algorithm  f 30 }  to  solve  the  unsteady  flame  equations.  They  used 
simplified  expressions  for  the  diffusive  transport  coefficients  vhich  were 
adjusted  to  give  good  laminar  flame  speed  prediction  for  methane-air  flames 
D1  xor.-Levis  has  more  recently  used  a  composite  flux  method  |3l!  to  stud 
a  variety  of  problems  like  the  kinetic  mechanism,  structure  and  propagaticr. 
of  flames  in  hydrogen— oxygen-nitrogen  flames  ;32j  and  flame  inhibition  by 
organic  halogen  compounds  1 3 2 i  •  He  discusses  the  ranges  of  applicability  cf 


he  oartia-  equilibrium  and  quasi-steady  assumptions  in  relation  t 
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tributicr.  cf  radical  populations  in  the  flames  [32]  • 

Warnatz  has  extensively  studied  both  freely  propagating  and  burner 


stabilized  flames  in  a  variety  of  premixed  gases  13^,  35,  3o,  3",  33.’ .  He 
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linearizes  the  chemical  reaction  t ernes  and  solves  the  time— dependent  equa¬ 
tions  implicitly  vith  assumed  initial  guesses  for  the  t erme nature  and  spec: 
orofiles.  He  also  uses  a  simplified  transport  model  [2~i  which  agrees  vel 
vith  the  ccmolete  multi cooocent  formulation  of  diffusion  and  thermal  coni' 


tier.,  .n  nos  recent  stu 
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he  concentration,  pressure  an:  tempera 


**  —  V  a  a  — o  'fo1 


ture  aeoer.aencfi  ox  trie  :  _a~e  veiocx. 


a  hydro gen -oxygen -nitrogen  mixtures 


concludes  that  at  the  present  state  of  knowledge  predictions  of  laminar, 
premixed  flame  propagation  should  be  as  reliable  as  experimental  results. 

Coffee  and  Heimerl,  using  a  method  of  lines  approach  [2?)  examined  di 
ferent  methods  of  approximating  miltispecies  transport  phenomena  in  models 
premixed,  laminar,  steady-state  flames.  They  concluded  coat  the  selection 
the  input  values  for  the  individual  species  transport  properties  is  more 
important  than  the  selection  of  the  approximation  method  !-C.. 

There  have  also  been  a  number  of  approaches  which  solve  a  steady  stat 
formulation  of  the  flame  equations.  These  include  the  works  of  Cixon-Levi 
[il]  ,  Wilde  \'-2\  and  Smooke  1^3) .  Che  advantages  and  the  difficulties  in 
using  steady  state  solution  procedures  have  been  discussed  by  Smooke 
These  methods  are  good  for  obtaining  burning  velocities  and  steady  state 
profiles,  but  cannot  provide  information  on  the  ignition  and  development  c 
flames. 

The  model  presented  in  this  report  selves  the  time -dependent  partial 
differential  equations  and  hence  can  describe  the  initiation,  propagation 
quenching  cf  laminar  flames.  It  uses  detailed  chemical  kinetic  models 
without  any  stead;.’  state  or  quasi— equilibrium  assumptions  and  sc  can  be  us 
to  stud;*  all  laminar  premixed  flames.  It  uses  an  asymptotic  coupling  net h 


in  conjunction  with  timestep  splitting  tc  couple  the  various  trysical  an: 
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i~ers,~ec.  “o  2sro  using  a.  cuaura.nins.-.-L^'  c cr/^rcsr.*  isn.icit  solution  c. 
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iteration  the  analytic  derivative  dA/d?  is  used 
ccmzutational  cell.  Thus 
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:i2). 

ATINC  assumes  that  pressure  and  density  are  constant  vithi.t  each  indi¬ 
vidual  finite-difference  cell  and  that  the  physics  is  evolving  slovly  enough 


for  full  communication  across  that  cell  to  have  occurred  in  a  times 


•vner. 


densities  in  adjacent  cells  are  different,  the  pressure 
two  cells  -ill  inn  art  =.  u—ier-su. 


gracter.t  betveer. 


right  anc  to  the  left  of  the  interface  betveen  the  tvc  ceils.  If  the  fluid 
vere  permitted  to  move  according  to  these  distinct  accelerations,  the  fluid 
•would  either  overlap  or  a  gap  -would  appear  at  the  interface  after  a  short 
•while.  To  prevent  this,  AT INC  matches  the  acceleration  calculated  from  the 
left  to  that  from  the  right.  This  is  described  in  more  detail  by  Boris  [l]. 
Problems  can  be  set  up  in  one  of  four  geometries:  cartesian  coordi¬ 
nates,  cylindrical  coordinates,  cylindrical  coordinates,  spherical  coordi¬ 
nates  and  a  variable  cover  series  coordinates  for  treating  one-dimensional 
nozzle  like  geometries.  Toe  boundary  conditions  are  controlled  by  specifying 
the  position  and  velocity  of  the  first  and  last  cell  interface. 
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subject  to  the  constraints 


,  =  2,  and  w  p  ,V_. 
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Ar.  iterative  algorithm,  DFI/JX,  enables  us  to  obtain  the  diffusion  velocitie 
vithcut  the  cost  of  performing  matrix  inversions  '2,3;  •  This  method  is  of 
2 O'!2';  and  is  vectorised.  Thus  it  is  substantially  faster  than  Z['A matrix 
inversions  when  four  or  more  species  are  involved. 

The  Zqs.  (l6,i"r)  are  conservatively  differenced  and  solved  explicit!;'. 
Solving  these  equations  requires  knowledge  of  the  binary  diffusion  coeffi- 

rn 

cients  2,.  ,  the  thermal  diffusion  ratios  K7  and  the  thermal 

o  *  <J 

conductivity  X.  Representations  of  these  coefficients  which  are  easily  use 
in  reactive  flow  models  and  ic  net  require  the  expensive  inversion  of  matri 
ces  are  chosen  and  these  are  discussed  in  a  subsequent  section  of  this 


'  1  — «oc*- 


ittins 


In  the  asymptotic  timestep  split  approach  '2j  ,  the  individual  terms  i: 
the  Tqs.  '1-1)  are  solved  independently  as  described  above  and  then  asynrtc 
ically  coupled  together.  Since  both  chemistry  and  sound  waves  are  usually 
stiff  in  deflagration  problems,  special  care  is  required  in  coupling  the 


hemical  heat  release  to  the  fluid  dynamics,  ur.  a  flame,  fluid  dynamic 


expansion  and  diffusive  transport  relieve  the  pressure  from  the  flame  regie 
as  fast  as  it  is  generated.  Thus  the  pressure  stays  effectively  constant. 
Small  pressure  fluctuations,  C(v_2/c2),  do  exist  and  are  just  large  enough  t 
drive  the  flows  which  reapportion  the  energy  released  by  chemistry-  or  trar.s 
ported  by  diffusive  processes. 
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sior.  and  thermal  conduction  are  converted  to 


cttve  ores  sure  cr.arees  as 
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oescrtoec  ce.ov.  oeicre  tne  flute.  cy ramie  dmestep,  tr.e  sue 
pressure  changes  is  added  to  the  old  pressure  distribution  to  obtain  a  rev 
pressure  distribution.  This  r.ev  pressure  distribution  is  used  to  modify  the 
ceil  er.trcpies.  The  r.ev  pressure  distribution  and  cell  entropies  are  used  in 
the  convection  transport  algorithm,  kZ I”C  to  obtain  the  convect  fluid  expan¬ 
sion  as  described  earlier. 

Tne  chemistry  step  should  be  taken  at  constant  pressure,  but  if  the 
profiles  change  slowly  per  cimestep,  it  nay  also  be  taken  at  constant  volume 
with  temperature  held  fixed.  On  completion  of  the  chemistry  integration,  tr.e 
heat  release  is  converted  to  an  effective  pressure  change  at  constant  volume 
because  the  cell  volume  has  been  held  fixed  during  the  chemical  kinetics 
calculation.  The  chemical  enery  released  car.  be  calculated  from  the  heats 
of  formation  and  the  changes  in  the  number  densities  of  the  various  species. 


Defining  an  enery  density  e  by: 


e  = 


-  foV2 


-  h  r. 


at  the  end  of  each  chemical  time  stem  we  have 


-t 
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Tr.e  changes  in  the  number  densities  are  calculated  by  solving  the  chemical 
rate  equations  as  described  in  an  earlier  section  of  this  report.  ’Jsir.g  f 
definition  of  enthalpy,  the  enery  density  s  is  written  as: 

r  s  r  v  c ' 
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afire ,  a  Sevton-Raphscn  iteration  technique  is  used  and  it  usually  converges 
in  one  or  tvc  iterations.  This  is  the  tenrerature  that  would  have  been 
attained  if  the  entire  energy  released  goes  into  heating  the  fluid.  Since 
the  chemical  energy  released  also  goes  into  fluid  expansion,  it  would  he 

use  this  as  the  temperature  for  the  next  timestep.  Therefore, 


the  energy  released  is  converted  to  an  effective  pressure  change  as 


IP  =  k_T  '  n  -  ? 


V  r 


hi' 


c^c 


where  T  is  the  temperature  calculated  from  So.  !22v  and  n,  are  the  new 


number  densities  of  the  various  species.  The  energy  transported  hy  thermal 
conduction  and  diffusion  are  also  converted  to  effective  pressure  changes  hy 
using  the  equation 

i?  *  hi  (y-d 

These  pressure  changes  are  used  as  energy  source  terms  for  the  fluid  dynamic 
t iciest et  as  described  above. 


.per,  zour.iary  ^ omit  ion 


In  the  study  of  unccnfir.ed  flames,  a  representation  for  the  open- 
bcur.dary  is  required  since  the  site  cf  the  computational  domain  is  limited. 
One  approach  is  to  allow  the  computational  tells  to  increase  in  sice  as  they 
’urther  from  the  flame.  Thus  the  comoutaticnal  domain  is  very  larce  ar. d 


set 


t.tere  is  no  corresponding  increase  in  computer  storage,  however,  the 


^  ^  ■*  r\ 


ihir.c  should  he  era  dual  to  limit  inaccuracies  which  arise  as  a  res-. 


tne  varying  cell  sice.  Furthermore,  adaptive  grinding  ’-“11  be  esse: 
this  approach,  since  as  the  flame  fror.t  moves  toward  the  large  cell: 
resolution  vill  have  to  decrease  around  the  flame. 

A  second  approach  which  we  have  used  to  stud;.’'  the  propagation  : 


.  w  ~  C.  _  n  _ 


fined  flames  takes  advantage  of  the  observation  that  the  effect 


o:  ar.  ooer. 


boundary  is  to  keep  the  pressure  essentially  constant.  Therefore, 
approach  the  movement  of  the  cell  boundary  that  is  on  the  open  end 
lated  by  allowing  the  pressure  in  each  cell  at  the  end  of  a  timeste 
adiabaticaily  to  the  oressure  before  the  timester.  In  oractice  thi 


cs 


ca_cu 
o  rel 
,s  cone 


by  calculating  the  changes  in  the  cell  volume  7„ ,  by 


I'd 
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■be 


P  =  ?(t)  +  A? 
new 


\2i', 


and  (A?)  is  the  change  in  the  presure  due  to  diffusive  transport,  energy 
released  in  chemical  reactions  and  any  energf  addition  to  the  system.  The 
changes  in  the  volume  of  the  cells  causes  the  location  of  the  open  boundary 
to  change.  The  location  of  the  open  boundary  and  the  fluid  velocity  in  the 
last  cell  (which  is  also  the  velocity  of  the  open  boundary)  are  used  as  open 
boundary  conditions  to  the  convective  transport  algorithm  ADIXC. 
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where  X?  is  tie  thermal  conductivity  of  gas  j  if  it  is  monatomic,  cr 
it  is  tie  thermal  conductivity  of  a  polyatomic  gas  with  internal  degrees 
freedom  "f  rccer." ,  and  M.  is  tie  atomic  mass  of  species  j.  one  thermal 

u 

conductivity,  X,  of  a  pure  polyatomic  gas  can  be  estimated  from  tie  non- 

v 

atomic  thermal  conductivities  X^  by  using  the  modified  Suker:  factor 
,  u, ,  as 

AJ=u,Xp  .2; 


E.  =  C.115  +  0.25^ 


:he  molecular  specific  heat  at  constant  pressure.  The  {X0, } 


may  ce  written  as 


, n  8. 322  x  io3  ,T  , 1/2 

-  *  77^' 
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u  0  %J 


[2  2^* 

vhere  1^”’  ,  a  collision  integral  normalised  to  its  rigid  sphere  value  is 

J  *J 

a  function  of  T*,  is  the  reduced  temperature  given  by 


The  terms  a,  and  s,  are  force  constants  in  tie  potential  function  which 

<J  v 

describes  the  interaction  between  two  molecules  of  the  same  species,  j. 


1.1.2  Ordinary  Diffusion  Coefficients 

We  assume  that  the  "ordinary"  (concentration)  diffusion  coefficients  for 


a  pair  of  species  (j,k)  ir.  a  multicocmonent  mixture  is 


:o  a  first  apprcxima- 


:hst  in  a  binary  mixture  of  species  J  and  k.  For  some  binary 


mixtures 


iilute  gases.  Mason  and  Marrero  [It;  have  gi 


ver.  a  semi-e: 
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;r  z he  calculation  of  the  thermal  conductivity,  the  bir.a ry  diffusi 


s  and  the  thermal  diffusion  ratios,  ve 
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potential  function  chosen.  For  the  results  presented  in  the  next  section,  ve 
have  used  the  ler.nard-Jones  {6-12}  potential  function  which  adequate!;'  de¬ 
scribes  the  interaction  of  non-tolar  molecules.  This  ootential  function  can 


■Titter  as 


p ' r }  =  It 


[if r2  -  iff 


where  e  is  the  detth  of  the  ootential  veil  (the  maximum  energy  of  attraction! 
and  o  is  the  collision  diameter  for  low  energy  collisions.  Tr.e  collision 
integrals  for  this  potential  function  have  been  calculated  and  tabulated  in 
many  independent  investigations  (e.g.  Ref.  [16],  [pc!).  Ve  have  used  the 
tabulations  of  Klein  et  al  [  p  3 ! •  There  is  some  uncertainty  in  the  values  of 
the  parameters  c  and  £  in  Zc.  (Ip).  Tr.e  parameters  ve  have  used  are  pre¬ 
sented  in  Table  II  and,  extent  for  <%-,  ou  n  and  euj  „  are  those 

“2^  **2^ 

given  by  Svehla  [pi].  Since  H20  is  a  polar  molecule  the  lennard-Jor.es  poten¬ 
tial  function  is  a  poor  approximation.  A  St ockaaye r-type  potential  function 
may  be  more  applicable  to  polar  molecules  [pi].  However  ve  have  assumed 
that  the  Ler.nard-Jones  potential  function  is  valid  for  H2C,  but  have  used  the 


values  giver,  by  Dixon-Levis  [33]  for  c-  «  an 


«  0  and  £—  -  (and  for  c-.} 

..2J  -2W 


which  he  has  found  to  be  satisfactory  in  extensive  investigations  of  the 
H 2-0 2  system. 

The  ordinary  diffusion  coefficients  have  beer,  calculated  using  Ic.  (33), 
where  the  values  of  A,,  and  have  been  taker,  from  the  data  giver. 


by  Mason  ar.d  Marrero  [1?[.  However,  r.c  such  data  is  available  for  many  pairs 
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Fig.  1  —Time  history  of  the  temperature  profile  in  a  2:1:10  mixture 
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the  cccmutation.s  with  a  good  guess  for  the  ter.per3.ture  and  species  profile 
=  re  ~icr.  ce'nind  —  fLane  frcnt#  closer  che  initial  profiles  ^re  t c  the 


:10a.,  steacy  state  :  -a~e  pro:  i_e ,  tr.e  scon* 


tr.e  stea: 


:pagating  flame. 


For  the  results  discussed  below,  the  model  was  set  up  in  Cartesian  gee 
etry  with  one  end  closed  'preventing  any  gas  flow  through  it;  and  the  other 
open  to  the  atmosphere.  At  the  closed  end,  a  first  approximation  to  the 
temperature  and  major  species  profile  was  set  up  as  initial  conditions. 

In  order  to  evaluate  the  imoortance  of  the  thermal  diffusion  in  the 
propagation  of  flames  we  have  studied  two  cases.  For  the  first  case  the 
effects  of  thermal  diffusion  have  net  been  considered,  that  is  the  source 


terms  due  tc  thermal  diffusion  in  uq.  t;  have  been  neglected. 


:ase  a, 


the  thermal  diffusion  ratio  H*.  cf  all  species  have  been  calculated  using 

w 

It.  (2dN»  However,  in  both  cases  the  heat  flux  due  to  concentration  gradi¬ 
ents  'the  Hu  four  effect)  has  not  been  considered  since  it  is  expect  tc  be 
negligible  even  when  thermal  diffusion  is  not  negligible  ’  1 0 ] . 


lase  1:  Thermal  Oiffusion  Teslected. 


Tr.e  calculations  were  initialized  as  described  above  and  soon 


a  crccara- 


ting  flame  developed  as  can  be  seen  in  Fig.  2,  where  the  temperature  profiles 


at  Ip  us  and  20  us  (from  the  sta 


**  -  - 


re  ccmutattons :  are  sr.cvr 
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required  since  the  flame  is  of  finite  thickness,  me  method  is  to  cheese  an 


a  criterion  for  the  -ccaticr.  of  tr.e  flame  is 


location  of  the  flame.  Tr.e  movement  cf  the  location  of 


arc  stare  md'.'es 
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Lty  :cmu:ed  using  IOC  grid  points  'cells!  is  ccriparec  tc  that  using 
^he  c c cmui at i on s  are  initialised  usinz  the  ur.iiDrml^" 
spaced  grid,  '-her.  the  number  of  grid  points  is  increased  the  spatial  res: 
tier,  increases  all  across  the  system  even  in  the  region  well  ahead  of  the 
flame  front. 


w'ase  2.  Effect  of  Thermal  Diffusion  on  fuming  velocity 

In  this  case  the  thermal  diffusion  source  terms  in  Eg.  (6)  were  included 
in  the  calculation  of  the  diffusion  velocities.  The  flame  velocity  and  the 
velocity  of  the  unturned  gases  were  determined  as  described  above  for  case  1. 
The  cime  history  of  these  velocities  have  been  compared  to  those  of  case  1  ir. 
Tig.  c.  The  effect  cf  thermal  diffusion  is  to  reduce  both  the  flame  velocity 
and  the  velocity  of  the  unburned  gases.  The  "steady"  flame  velocity  (corre¬ 
sponding  to  20  us  from  the  start  of  the  calculations)  decreases  from  1C  m/s 
to  about  9*65  m/s  while  the  corresponding  velocity  of  the  unburned  gases 
decreases  from  7. S3  c/s  to  about  "T. 65  m/s.  The  net  effect  is  to  lower  the 
burning  velocity  of  the  stoichiometric  hydrogen-air  mixture  to  about  2  m/s. 
The  reasons  for  these  effects  are  examined  below. 

The  source  terms  considered  in  Ec.  ("0  for  the  calculation  of  diffusion 
velocities  are: 

(a)  due  to  concentration  gradients  (ordinary  diffusion) 


S  =  V(-*0 
0D  ' 


'us) 


(b)  due  to  temperature  gradients  (thermal  diffusion) 


c  —  y-  *  * 
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The  scatial  variations  of  these  two  terms  have  been  calculated  for  the  steady 


TIME  (ms) 

.  8  —  Effect  of  thermal  diffusion  on  the  time  history  of  the  velocity  of 
(a)  the  flame  and  (b)  the  fluid  at  the  open  boundary 


propagating  flame  ccresponding  to  2C  us  in  Fig.  " .  The  temperature  profile 
across  the  system  am  2C  us  is  given  in  Fig.  9*  For  convenience,  the  abscissa 
is  given  in  cell  numbers.  In  Fig.  5,  the  spatial  variation  of  the  source 
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r.  compared  for  the 

species  C2  and  S2.  We 

see  that 

for 

C, ,  the  term 
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,s  opposite  in 

sign  to  5^r.  A  negative 

value  for 

Dm-' 

implies  that 

the 

oxygen  mol ecu! 

es  are  diffusing  from  left 

to  right ,  that  is,  from  the  hotter  region  to  the  cooler  region.  Since  the 
sum  of  the  two  source  terms  determines  the  diffusion  velocity,  the  effect  of 
thermal  diffusion  is  to  slow  the  diffusion  of  oxygen  molecules  into  the 
higher  temperature  region  of  the  flame.  The  effect  of  this  on  the  burning 
velocity  is  not  obvious  since  it  depends  on  the  details  of  the  reaction 
mechanism  and  the  distribution  of  other  radicals.  For  the  species  X,,  the 
term  3mn  is  comparable  to  the  term  and  the  effect  of  thermal 
diffusion  is  to  enhance  the  movement  of  N2  from  the  hotter  region  of  the 
flame  to  the  cooler  region. 

In  Fig.  1C,  we  compare  the  terms  for  the  species  rl2  and  K.  For 
H2,  the  effect  of  thermal  diffusion  is  similar  to  that  of  ordinary  diffusion 
and  it  enhances  the  movement  of  H2  from  the  cooler  region  to  the  hotter 
region.  The  term  3gr  for  the  species  H  is  particularly  interesting.  It 
is  observed  that  the  K  atoms  in  the  region  up  to  cell  number  To  (i.e.  T  > 
luCO  K  have  a  tendency  tc  diffuse  into  the  hotter  region  while  those  after 
cell  number  75  (T  <  14(30  K)  diffuse  into  the  cooler  region.  The  diffusion  of 
H  atoms  into  the  cooler  region  is  very  important  in  the  propagation  of  the 
flame.  The  effect  of  thermal  diffusion  is  to  delay  and  decrease  the  diffu¬ 
sion  of  H  atoms  into  the  cooler  region.  Since  the  H  atoms  are  produced  pri¬ 
marily  in  the  high  temperature  region  (in  fact  here  the  concentration  of  H 
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Fig.  11  —  Comparison  of  the  spatial  variation  of  the  thermal  diffusion  source  term  for  the 
species  H  and  0.  The  temperature  profile  across  the  system  has  also  been  depicted. 
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This  report  described  a  one-ciraensional,  tine  -dependent ,  lagrangian 
odel  for  predicting  the  initiation,  propagation  and  quenching  of  laminar 
lanes  in  pre-mixed  gases.  A  number  of  new  approaches  and  algorithms  used 
he  model  have  been  discussed.  An  asymptotic  coupling  method  was  used  in 
or. Junction  with  timestep  splitting  to  couple  the  various  physical  and  che 
al  processes.  This  approach  allowed  the  use  of  entirely  different  algo- 
ithms  for  the  processes  represented  by  different  terms  in  the  ccr.servatic 
quaticns.  This  made  the  computing  procedure  very  efficient  and  i.cexpen- 
ive. 

The  result  from  several  calculations  using  the  model  have  been  pre- 
er.ted.  The  first  is  a  flame  initiation  and  minimum  ignition  energv  stud:' 

hj'ircgen,  oxygen,  nitrogen  mixture.  This  used  the  time -dependent  proper 
f  the  model  and  the  obervations  were  in  qualitative  agreement  with  experi 
e.ntal  data.  Quantitative  comparisons  were  not  possible  since  the  compos i 
ion  of  the  mixture  and  the  time  for  ener®"  deposition  were  different.  Th 
econd  problem  discussed  was  an  estimation  cf  the  steady  state  burning  vel 
ty  of  a  stoichiometric  hydrogen-air  mixture.  The  burning  velocity  premie 
y  the  model  was  on  the  lower  end  cf  the  observed  experimental  range.  We 
cte,  however,  that  there  are  several  factors  which  would  tend  to  increase 
his  estimated  value.  First,  the  calculated  burning  velocity  would  be  lar 
f  the  velocity  of  the  unburned  gas  used  in  -his  calculation  was  closer  to 
he  flame  front  instead  of  at  the  c pen  boundary.  It  is  not  clear  where  it 
'could  be  evaluated  in  a  time-dependent  calculation.  Second,  the  ur.certai 
ies  in  the  values  cf  some  of  the  diffusion  coefficients  are  about  27T.  7 


id 


burning  velocity  was  found  to  be  very  sensitive  to  certain  diffusion 
coefficients.  Because  of  the  uncertainty  in  the  experimental  values  of  tee 
diffusion  coefficients,  there  nay  net  be  any  advantage  in  using  the  empirical 
formula  over  the  expression  from  kinetic  theory. 

The  model  also  oredicts  that  the  effect  of  thermal  diffusion  is  to  lover 


the  burning  velocity  by  abou 
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ir.is  observation  is  m  agreement  vitr.  tr.at 


of  Dixon-Levis  '22;  vho  found  that  neglecting  thermal  diffusion  of  hydrogen 
atoms  increases  the  burning  velocity  by  $  or  6'.  For  the  hydrogen  atom,  the 
thermal  diffusion  source  term  vas  found  to  be  small  compared  to  the  source 


;erm  due  tc  ordinary  (concentration)  diffusion.  3ut  because 
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chemistry  occurring  in  flames,  a  small  difference  in  the  hydrogen  source  term 
causes  a  significant  difference  in  the  burning  velocity.  This  again  empha- 
sizes  the  need  for  more  accurate  estimates  of  the  diffusion  coefficients. 

In  conclusion,  the  model  can  be  used  confidently  to  predict  the  igni¬ 
tion,  propagation  and  quenching  of  laminar  flames  in  pre-mixed  gases  provided 
there  is  satisfactory  data  on  the  chemical  kinetics  and  the  diffusive  trans¬ 
port  coefficients. 


ACKNOWLEDGMENTS 

The  authors  would  like  to  acknowledge  the  help  and  advice  of  T.E.  Young, 
J.M.  Picor.e,  ar.d  Homer  Carhart.  This  work  is  sponsored  by  the  Office  of 
Naval  Research  and  the  Naval  Material  Command. 


-TTr"3r”'"rc 


=cr.s, 


Imclicit  laeranriar.  Hvirodyr.amics  Tods,  NR.l 


Memorandum  Pep  or"  -C22,  Naval  Research  laboratory ,  Washington ,  1C,  Ip 
Bran ,  Z.3.,  arc  lords,  J.P. ,  ?rcg.  Energy  B  oncost  lor.  5cier.ce,  ~ : l-~2 


’ones,  W.V,  ,  and  Boris,  J.P. 


JT.Sin*  ~ 


Yeung,  T.R.  ,  and  Boris,  J.P.  ,  J.  Rhys.  Chen,  dl,  2-2^-2~2~  [ip’"’,'. 
Young,  T.R.,  CHZMZQ  -  A  Subroutine  for  Solving  Stiff  Ordinary 
Differential  Equations,  M?.l  Memo  ran  dus  Report  -091,  Naval  Research 
laboratory,  Washington,  DC,  I960. 

Bran,  E.A. ,  and  Boris,  J.P. ,  presented  at  the  "th  International 
Colloquium  on  Das  Dynasties  of  Explosions  and  Reactive  System, 

M  _ 

Goettingen,  19~9»  appears  in  rrog.  in  Astronautics  and 
Aeronautics,  *'6:1;1-1"1  (1951). 

Eailasanath,  K.  ,  Bran,  E.3.,  and  Boris  J.P.,  A  Theoretical  Study  of  z 
Ignition  of  Pre-Mlxed  Gases,  lyfi2  (to  appear  in  Combust.  Plane). 
Irdritz,  D. ,  Boris,  J. ,  Carhart,  H. ,  Iran,  E. ,  Sheinson,  R. ,  Williams 
F. ,  and  Young,  T. ,  Computation  of  Hydrogen-Oxygen  Flammability  limits 
submitted  to  Comb.  Flame. 

Eailasanath,  ,  Oran,  E.3.,  Boris,  J.P.,  and  Young,  T.R.,  Bime- 
Bependent  Simulation  of  Flames  in  Hyircger-Cxyren-Nitroger.  Mixtures, 
appear  in  the  Proceedings  of  the  GAMM  Workshop  on  Numerical  Methods  i 
laminar  Flame  Propagation,  Vieveg-Verlag,  w.  Germany,  IJcl. 

Williams,  F.A.,  Combustion  Theory.  Addison  Wesley,  Reading,  MA ,  1965, 


-  -j*  *- 


V-  - 


51.  Stull,  1.?..,  and  Prophet,  H. ,  DANAF  Thermechemical  Tables ,  National 
Standard  Reference  lata  Series,  'National  Bureau  of  Standards,  1c.  I ~ , 

2nd  Ed. ,  Gaithersburg,  ME,  1?~1. 

“2.  Oran,  E.S.,  Young,  T.R.,  Boris,  J.P.,  and  Cohen,  A.,  Weak  anc  Strong 

Iscition-I.  Numerical  Simulations  of  Shock  Tube  Ercoeriments .  NRL  Memo¬ 
randum  Report  1664,  Naval  Research  Laboratory,  Washington,  DC,  1981. 
(Also  to  appear  in  Combust.  Flame). 

53.  Klein,  M. ,  Hanley,  H.J.M.,  Smith,  F.J. ,  and  Holland,  ?. ,  Tables  of 
Collision  Integrals  and  Second  Virial  Coefficients  for  the  ln,6,3) 
Ir.temolscular  Potential  Function,  National  Standard  Reference  Data 
Series,  National  Bureau  of  Standards,  No.  L"7,  Gaithersburg,  ME,  19"-. 

-.  Svehla,  R.A. ,  Estimated  Viscosities  and  Thermal  Conductivities  of  Gases 
at  High  Temperatures ,  Technical  Report  No.  R-L32,  NASA,  Washington,  DC, 
1962. 

35.  Levis,  3.,  and  von  Elbe,  G. ,  Combustion,  Flames  and  Explosions  of  Gases, 
.Academic  Press,  NY,  19cl,  pp.  323-36'7. 

56.  Baulch,  D.L.,  Drysdale,  D.C.,  Horne,  D.C.,  and  Lloyd,  A.C.,  Evaluated 
Kinetic  Data  for  Hirh  Temperature  Reactions,  Yoi.  1,  Buttsrvcrths , 
London,  1?"2. 

3".  Sampson,  R.F.,  and  Garvin,  D. ,  Chemical  Kinetic  and  Photochemical  Data 
for  Modelling  of  .Atmospheric  Chemistry,  NBS  Technical  Note  366,  'National 


Bureau  cf  Standards,  Washington,  DC  '19"5). 


3c.  Cohen,  N. ,  and  Westoerg,  K.R. ,  Data  Sheets,  The  Aerospace  Dorpc ration 


Box  ?2?5",  les  .Angeles,  DA  (l?"?;, 


